home *** CD-ROM | disk | FTP | other *** search
/ Turnbull China Bikeride / Turnbull China Bikeride - Disc 1.iso / ARGONET / PD / MATHS / RLAB / RLAB125.ZIP / !RLaB / toolbox / quadr < prev    next >
Text File  |  1995-02-19  |  5KB  |  180 lines
  1. % QUADR Numerical evaluation of an integral.
  2. %    Q = QUADS(F,X1,X2) approximates the integral of F(X)
  3. %    from X1 to X2 with relative error  1e-4.  'F' is a string
  4. %    containing the name of the function such as 'sin' or expression
  5. %    depending on x (such as '1/(x+1)^2'). Function or expression F
  6. %    must return a vector of output values of the same size as a
  7. %    vector of input values.
  8. %
  9. %    Q = QUADS(F,X1,X2,TOL) integrates to a relative error of TOL.
  10. %    Q = QUADS(F,X1,X2,TOL,TRACE) also traces the function
  11. %    evaluations with a point plot (with different colors for each
  12. %    recursion level).
  13. %    Q = QUADS(F,X1,X2,TOL,TRACE,P1,P2,...) allows parameters P1, P2
  14. %    to be passed directly to a function F in the form F(X,P1,P2,...).
  15. %    When F stands for expression, its parameter dependence must
  16. %    contain 'p1' ... explicitly, such as '(x-p1)/(x+p2)^p3'.
  17. %    [Q,X,Y] = QUADS(F,...) also returns vector of points X of
  18. %    function evaluation and corresponding function values Y.
  19.  
  20. %  Kirill K. Pankratov,  kirill@plume.mit.edu
  21. %  01/30/95
  22.  
  23. require lininsrt
  24.  
  25. quadr = function (fun, x1, x2, tol, trace, p1, p2, p3, p4, p5)
  26. {
  27.   local (x1, x2, tol, trace, p1, p2, p3, p4, p5)
  28.  
  29.   % Defaults and parameters ..................................
  30.   tol_dflt = 1e-4; % Default tolerance for integral estimation
  31.   trace_dflt = 0;  % No plotting by default
  32.   n_max = 7;       % Max. number of points inserted at one
  33.                    % iteration between 2 old points
  34.   max_lev = 8;     % Max. number of iterations (recursion level)
  35.   pow = 1/4;       % Power index for estimation of number of
  36.                    % inserted points
  37.  
  38.   % Coefficients ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  39.   % Function interpolation coefficients ..........
  40.   % (Lagrange polynomials interp. from 4 basis pts to 3 pts.
  41.   ci =  [-1/6, 2/3,  0,  2/3, -1/6];   % [1 2 4 5] to 3
  42.  
  43.   % Integral approximation coefficients ..........
  44.   cq = [0.31111111111111;
  45.         1.42222222222222;
  46.         0.53333333333334;
  47.         1.42222222222222;
  48.         0.31111111111111];
  49.  
  50.  
  51.   % Handle input .............................................
  52.   if (nargs < 5) { trace = trace_dflt; }
  53.   if (nargs < 4) { tol = tol_dflt; }
  54.   if (trace == []) { trace = trace_dflt; }
  55.   if (tol == [] || all([tol[:]; -1]) <= 0) { tol = tol_dflt; }
  56.  
  57.   % Set up function call .....................................
  58.   if (class (fun) == "function")
  59.   {
  60.     call = "fun ( x";
  61.     for (n in 1:nargs-5)
  62.     {
  63.       call = call + ", p" + num2str(n);
  64.     }
  65.     call = "f = " + call + ");";
  66.  
  67.   else               % .. Expression
  68.  
  69.     call = "f = " + fun + ";";
  70.   }
  71.  
  72.   % Initial layout ..............................
  73.   g = sqrt(3)/3;  % Initial layout section ratio (just some
  74.                   % number close but not equal to 1/2)
  75.   span = x2 - x1;
  76.   x_o = [x1, x1+[1-g, 1+g/117, 1+g]*span/2, x2];
  77.  
  78.   % Auxillary ...................................
  79.   ci5 = ci[5];
  80.   ci = ci[1:4];
  81.   cq5 = cq[5];
  82.   cq = cq[1:4]';
  83.   tol_f = tol/sqrt(span)*16;
  84.   o4 = ones(4,1);
  85.  
  86.   % Initialize output and control parameters ....
  87.   f_o = [];
  88.   n_ins = 3;
  89.   n_i = 1;
  90.   last = 1;
  91.   iter = 0;
  92.   is_more = 1;
  93.  
  94.   % Begin recursively refine approximations ..................
  95.   while (is_more)  % While not enough accuracy ``````````````````0
  96.   {
  97.     iter = iter + 1;
  98.  
  99.     % Insert new points ..........................
  100.     F = n_ins[o4;];
  101.     </ mask ; x_o /> = lininsrt(x_o, F);
  102.     if (iter == 1) { mask = ones(size(mask)); }
  103.  
  104.     % Expand "last" vector
  105.     last = last[o4;1:n_i];
  106.     last = [last[:]', last[n_i]];
  107.  
  108.     % Evaluate function at new points ............
  109.     x = x_o[find (mask)];
  110.     eval(call);
  111.  
  112.     % Insert new points in the output and control vectors
  113.     n_ins = find(!mask);
  114.     misfit = find(mask);
  115.     if (iter > 1)
  116.     {
  117.       f_o[n_ins] = f_o;
  118.       last[n_ins] = last;
  119.     }
  120.     f_o [misfit] = f;
  121.     last[misfit] = ones(size(misfit));
  122.  
  123.     % Reshape function values vector ............
  124.     l_o = length(x_o);
  125.     n_i = (l_o - 1)/4;
  126.     F = reshape(f_o[1:l_o-1], 4, n_i);
  127.  
  128.     % Calculate difference between actual 
  129.     % and interpolated values ...................
  130.     misfit = ci*F+ci5*f_o[5:l_o:4];  % Interpolation
  131.     misfit = abs(misfit-F[3;]);      % Difference
  132.  
  133.     % Calculate how many points need to be inserted
  134.     scale = mean(abs(f_o));
  135.     n_ins = (misfit/(tol_f*scale)).^pow-1;
  136.  
  137.     % Smooth (moving average)
  138.     if (iter > 2)
  139.     {
  140.       misfit = cumsum(n_ins);
  141.       n_ins[3:n_i-2] = (misfit[5:n_i] - misfit[1:n_i-4])/4;
  142.       n_ins[[1, 2, n_i-1, n_i]] = n_ins[3, 3, n_i-2, n_i-2];
  143.     }
  144.     
  145.     % Make it integer and upper bounded
  146.     n_ins = min(ceil(n_ins), n_max*ones(size(n_ins)));
  147.  
  148.     % If previous iteration was accurate
  149.     % enough, do not insert new points .........
  150.     misfit = n_ins > 1;
  151.     n_ins = max(n_ins, last[5:l_o:4]);
  152.     last = misfit;
  153.  
  154.     is_more = any(n_ins);   % If needs more points inserted
  155.  
  156.     if (trace)    % Plot new points ..............
  157.     {
  158.       plot([ x, f ]);
  159.     }
  160.  
  161.     % If still not enough accuracy .............
  162.     if (iter > max_lev)
  163.     {
  164.       disp(" Recursion level limit reached. Singularity likely.");
  165.       is_more = 0;
  166.     }
  167.  
  168.   }  
  169.   % End while (recursion) ''''''''''''''''''''''''''''''''0
  170.  
  171.   if (trace) { "hold off" }
  172.  
  173.   % Now calculate the integral itself .........
  174.   last = 5:l_o:4;
  175.   misfit = cq*F+cq5*f_o[last];
  176.   Q = (x_o[last] - x_o[last - 1])*misfit';
  177.  
  178.   return << Q=Q; x_o=x_o; f_o=f_o >>;
  179. };
  180.